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ABSTRACT 

We have used data from the Sloan Digital Sky Survey (SDSS) Data Release 5 to explore the overall structure 
and substructure of the stellar halo of the Milky Way using ^ 4 million color-selected main sequence turn-off stars 
with 0.2 < g-r < 0.4 and 18.5 < r < 22.5. We fit oblate and triaxial broken power-law models to the data, and 
found a 'best-fit' oblateness of the stellar halo 0.5 < c/a < 0.8, and halo stellar masses between Galactocentric 
radii of 1 and 40kpc of 3.7 ± 1 .2 x 1O**M0. The density profile of the stellar halo is approximately p cx r"", where 
— 2 > a > -4. Yet, we found that all smooth and symmetric models were very poor fits to the distribution of stellar 
halo stars because the data exhibit a great deal of spatial substructure. We quantified deviations from a smooth 
oblate/triaxial model using the RMS of the data around the model profile on scales > 100 pc, after accounting for 
the (known) contribution of Poisson uncertainties. Within the DR5 area of the SDSS, the fractional RMS deviation 
(j/total of the actual stellar distribution from any smooth, parameterized halo model is > 40%: hence, the stellar 
halo is highly structured. We compared the observations with simulations of galactic stellar halos formed entirely 
from the accretion of satellites in a cosmological context by analyzing the simulations in the same way as the 
SDSS data. While the masses, overall profiles, and degree of substructure in the simulated stellar halos show 
considerable scatter, the properties and degree of substructure in the Milky Way's halo match well the properties 
of a 'typical' stellar halo built exclusively out of the debris from disrupted satellite galaxies. Our results therefore 
point towards a picture in which an important fraction of the stellar halo of the Milky Way has been accreted from 
satellite galaxies. 

Subject headings: Galaxy: halo — Galaxy: formation — Galaxy: evolution — galaxies: halo — Galaxy: 
structure — Galaxy: general 



1. INTRODUCTION 

The stellar halo of the Milky Way has a number of distinctive 
characteristics which make it a key probe of galaxy formation 
processes. Milky Way halo stars have low metallicity, alpha 
element enhancement, a high degree of support from random 
motions, and a roughly r"^ power law distribution in an oblate 
halo (Eggen, Lvnden-Bell & Sandage 1962; Chiba & BeersI 
I2OOO; Yannv et al. 2000; Larsen & Humphreys 2003; Lemon 
et al. 2004; Newberg & Yanny 2005; Juric et al. 2007). The low 
metallicities and alpha element enhancements suggest that the 
stars formed relatively early in the history of the Universe. Yet, 
there has been disagreement about where these stars formed: 
did they form in situ in the early phases of the collapse of the 
Milky Way (e.g., Egg en. Lvnden-Bell & S andage 1962). or did 
the stars form outside the Milky Way in satellite g alaxies only to 
be acc reted by the Milky W ay at a la ter date (e.g..l Searle & Zing 
[19781; iMajewski, Munn, & Hawleylll 996; Bullock, Kratsov, & 
Weinberg 2001; Bullock & Johnstonli2005i; Abadi, Navarro, & 
Steinmetz .2006.) ? 



A key discrimin ant between thes e pictures is the structure 
of the stellar halo ( Maiewskilll993l) . In situ formation would 
predict relatively little substructure, as the formation epoch 
was many dynamical times ago. In contrast, current models 
of galaxy formation in a hierarchical context predict that the 
vast majority of stellar h alo stars should be accreted from dis- 
rupted satellite galaxies (Johnston' 1998; Bullock, Kratsov, & 
Weinberg 2001; Bullock & Johnston 2005; Moore et al. 20QS 
I Abadi. Navarro. & Steinmetz! 120061) . The accumulated debris 
from ancient accretion episodes rapidly disperses in real space 
(although in phase space, some informat ion about initial condi- 
tions persists; e.g..' Helmi & Whit3ll999 ). forming a relatively 
smooth stellar halo. The debris from accretions in the last few 
Gyr can remain in relatively distinct structures. Simulations 
predict quite a wide range in 'lumpiness' of stellar halos, with 
a general expectation o f a significant amount of reco gnizable 
halo substructure (B ullock. Kratsov. & Weinberg|200U Bullock 
& Johnston 2005 ). 

Consequently, a number of groups have searched for sub- 
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structure in the Milky Way's stellar halo, identifying at least 
3 large-scale features — tidal tails from the disruption of the 
Sagittarius dwarf galaxy, the Low-L atitude stre am, and the 
Virgo overdensity (llbata. Gilmore. & Irwin 1995; Yanny et alj 
2000; Ivezic et al. 2000; Newberg et al. 2002; Maiewski et al. 
2003fc Jannv et al. 2003; Ibata et al. 2003; Juric et a L2007i; D uf- 
fau et al. 2006; Belokurov et al. 2006a; Newberg et al.ll2007l al- 
though see Momany et al. 2006 for a discussion of possible dis- 
rupted disk origin of much of the Low-Latitude stream) — and 
a ho st of tidal tails from globular clusters (e.g.. Odenkirchen 
et al. '2003'; 'Grillmair & Johnsori'2006'), dwarf galaxies (e.g., 
llrwi n & Hatzidimitriou 1995; Martmez-Delgado et al. 2001), 
and of unknown origin (e.g ., Belokurov et al. 2006b; G rillmairi 
l2006al ; lGrillmair & Dionatos.2006.;.Belokurov et al. 200 f). Fur- 
thermore, substruct ure has been observed in the st ellar halos of 
other galaxies (e.g..'S hang eTalll 19981; llbata et al.ll2001i) . Thus, 
it is clear that accretion of stars from satellite galaxies is a con- 
tributor to the stellar halos of galaxies. 

Yet, it remains unclear whether accretion is the dominant 
mechanism for halo build-up. A key observable is the fraction 
of stars in substructure (or a quantitative measure of the degree 
of substructure): if much of the halo mass is held in substruc- 
tures, this argues for an accretion origin; if instead a tiny frac- 
tion of halo stars is held in substructures, this places (very) tight 
constraints on any recent accretion scenario. However, it is not 
clear how best to address this question. One possible approach 
is to define 'overdense' areas of the halo by hand or algorithmic 
means, and to fit the rest with a smooth hal o component; the re- 
maind er would be in 'overdensities' (e.g.. iNewberg & Yannvl 
120051) . Here, we take a different approach. Since one does 
not know a priori which stars should be 'smooth halo' stars 
and which are in 'overdensities', we treat all halo stars equally, 
fit a smooth model, and examine the RMS of the data around 
that smooth model (accounting for the contribution to the RMS 
from counting statistics). In this way, we obtain a quantitative 
measure of the degree of halo structure on > lOOpc scales with- 
out having to make uncomfortable decisions about which stars 
should be fit with a smooth component and which should be 
included in overdensities. 

In this paper, we apply this technique to explore the struc- 
ture of the stellar halo of the Milky Way, and place constraints 
on the fraction of stars in stellar halo under- or over-densities 
using imaging data from the Fift h Data Release ( DR5) of the 
Sloan Digital Skv Survev (SDSS: 'York et al. 2000: Adelman- 
McCarthy et al. 2007). Under the assumption that the bulk of 
the stellar population in the stellar halo is relatively metal-poor 
and old, we isolate a sample dominated by halo main sequence 
turn-off stars and explore the distribution of halo stars as a func- 
tion of Galactic latitude, longitude and distance from the Sun 
(^. In |3] we generate a grid of smooth halo models to com- 
pare with the data, and in ^we constrain the 'best-fit' smooth 
stellar halo parameters and quantify the fraction of halo stars in 
stellar halo under- or over-densities. We compare the observa- 
tions with models of stellar halo formation in a cosmological 
context in ^ and present a summary in ^ 

2. DATA 

SDSS is an imaging and spectroscopic survey that has mapped 
~ 1 /4 of the sky. Imaging data are produced si multaneously in 
five photometric bands, namely u, g, r, i, and z (Fukugita et al.l 
IT996I: iGunn et alJ[T998l: iHogg et alJl200ll : iGunn et al.ll2006E 
The data are processed through pipelines to measure photomet- 
ric and astrometric properties (.Lupton. Gunn. & Szalav.. 19991 : 



IStoughton et alJ2002tlSmith et al.ll2002tlPier et alJl2003l: Ivezic 
et al. 120041 : iTucker et alj|2006t) and to select targets for spectro- 
scopic follow-up. DR5 covers ~ 8000 square degrees around 
the Galactic North Pole, together with 3 strips in the Galactic 
southern hemisphere. We use the catalog of objects classified 
as stars with artifacts removedQ, together with magnitude limits 
r < 23.5 and g < 24.5. Photo metric uncertainties as a function 
of magnitude are discussed in ISesar et"a D (|2007). We choose 
to analyze only the largely contiguous ^ 8000 square degree 
area around the Galactic North Pole in this work, giving a to- 
tal sample of ^ 5 x 10^ stars, of which ^ 3.6 x 10'' stars meet 
the selection criteria we apply later. In what follows, we use 
Galactic extinction c orrected magnitudes and colors, following 
ISchlegel et all [l9m : such a correction is appropriate for the 
stars of interest in this paper owing to their large heliocentric 
distances ZJheiiocentric > 8kpc. 

2. 1 . Color-magnitude diagrams: an introduction 

To help get one's bearings, it is instructive to examine some 
color-magnitude diagrams (CMD) derived from these data (Fig. 
[TJ. The color-magnitude diagram of all stars with b > 30° is 
shown in the top left panel, where the grey levels show the log- 
arithm of the number of stars in that bin per square degree from 
lO"-' stars/deg^ to 7.1 stars/deg^; such a scaled CMD is fre- 
quently called a Hess diagram. To help interpret this Hess dia- 
gram, we show two additional Hess diagrams for two globular 
clusters covered by these data: Pal 5 and NGC 5024 (in what 
follows, distances and metallici ties for these and all other glob- 
ular clusters are adopted from iHanisI 1 1 9961) . The top middle 
panel of Fig. [T] shows a Hess diagram for stars in the globular 
cluster Pal 5 (a circle of radius 0.°5 around the position / = 0.°85 
and b = 45 .°9). The grey levels show: 

(A^onA;^-A^offA;^)/^offA-^, (1) 
where denotes the number of stars in the field of interest (de- 
noted by the subscript 'on') and a control field 'off, and A is 
the area of that field. In this case the control field is nearby: 
a circle of radius 4° around the position / = 6° and b = 46°. 
One can clearly see the main sequence turn off with g-r ^ 0.2 
and r ^ 21, with the lower main sequence extending redwards 
towards fainter magnitudes and the subgiant branch extending 
redwards towards brighter magnitudes. In the top right panel 
we show a similar Hess diagram for NGC 5024; because this 
cluster is rather brighter than Palomar 5 the CMD is better pop- 
ulated and shows a more prominent red giant branch (extend- 
ing towards brighter magnitudes with g-r ^ 0.5) and horizonal 
branch (with g-r <0 and r ^ 17). 

There are a few points to note about Fig. [T] Firstly, for 
old populations such as those in globular clusters it is clear 
that the color of the main sequence turn-off (MSTO) is a 
metallicity indicator (this point is discuss ed in more detail for 
SDSS isochrones in iGirardi et al.ll2004 . Comparing Pal 5 
([Fe/H] ~ -1.4, (g-r)MSTO ~ 0.3) with NGC 5024 ([Fe/H] - -2.1, 
{g-r)MSTO ^ 0.15), one can see that old very metal-poor popu- 
lations ([Fe/H] < -2) have bluer main sequence turn-offs com- 
pared to less metal-poor populations ([Fe/H] ~ -1.5). Sec- 
ond, MSTO stars are a reasonably good distance indicator, al- 
beit with significant scatter In Fig. |2] we show the absolute 
magnitude distribution of all stars with 0.2 < g-r < 0.4 in 
Pal 5 (solid fine: distance= 22.6 kpc), NGC 5024 (dashed line: 
distance= 18. 4 kpc) and a third globular cluster NGC 5053 (dot- 
ted line: [Fe/H] ^ -2.3, distance= 16.2kpc). The mean r-band 

'See |http : //cas ■ sdss ■ org/astro/en/help/docs/realquery . asp#f 1 




Fig. 1 . — Hess diagrams in terms of g — r color and r-band magnitude derived from tlie SDSS data. In these Hess diagrams, we sliow for completeness the data to 
the very faintest limits r > 23, where the S/N is low and there is significant contamination by misclassified galaxies, spurious detections, etc. These diagrams show 
in general two plumes in the stellar density distribution that reflect main sequence turn-off stars with g-r ^ 0.3 and intrinsically faint and low-mass disk dwarf 
stars with g—r ~ 1.4. We limit our analysis to 18.5 < r < 22.5 in this paper for the main sequence turn-off dominated color bin 0.2 < g—r < 0.4, in the area where 
the data quality is still excellent. Top left: The density of stars per square degree per color interval per magnitude for h > 30°, scaled logarithmically. This Hess 
diagram contains 4 X lO' stars. Top middle: The Hess diagram for the (sparsely-populated) globular cluster Pal 5 (within a circle of radius 0?5). Top right: The 
Hess diagram for the globular cluster NGC 5024. Bottom left: A difference Hess diagram (following Eqn.[T) differencing two lines of sight (l,h) = (300,70) and 
{l,b) = (60,70). The grey scales saturate at ±100%. In an axisymmetric halo, this difference should equal zero within the shot noise: it clearly does not. Bottom 
middle: A difference Hess diagram differencing two lines of sight {l,h) = (44,40) and (l,b) = (15,45). The grey scales saturate at ±50%. This Hess diagram should 
be close to, but not exactly equal to, zero. Bottom right: A difference Hess diagram differencing two lines of sight (l,b) = (167,35) and (/,&) = (193,35). The grey 
scales saturate at ±50%. Again, in a symmetric halo, this difference should equal zero. 



absolute magnitudes of the distributions are (4.3,4.7,5.0) re- 
spectively, and all distributions individually have RMS ~ 0.9 
mag. Thus, modulo a metallicity-dependent < 0.5 systematic 
uncertainty, the MSTO is a good distance indicator with ~ 0.9 
mag scatter. 

Examining the top left panel of Fig.lT] in the light of the glob- 
ular cluster CMDs, it is possible to interpret some of the fea- 
tures of the b > 30° Hess diagram. At all distances, the MSTO 
is visible as a clearly-defined as a sharp 'blue edge' to the dis- 
tribution, indicating to first order that the stars in the galactic 
disk at large scale heights and in the stellar halo are dominated 
by a metal-poor old population with ages not that dissimilar to 
those of globular clusters; this is the assumption that we will 
adopt in the remainder of this paper At g- r < 0.5, one sees 
the MSTO for stars in the stellar disk at >kpc scale heights (at 
r < 18; often the disk at such scale heights is referred to as the 
thick disk) and in the stellar halo (at r > 18). One can see a 
'kink' in the MSTO at r ^ 18, as highlighted by the contours; 
we interpret this as signifying a metallicity difference between 
the disk at ~ kpc scale heights and stellar halo (this transition 



is also ver y prominent in Fig. 4 of iLemon et al. 1 12004 ' and in 
IChen et al . 2001, who interpret this CMD feature in the same 
way). The MSTO in the stellar halo has g-r ^ 0.3, similar 
to that of Pal 5 ([Fe/H] ---1.4) and 0.15 mag redder than 
those of NGC 5024 and NGC 5053 with ([Fe/H] < -2). This 
suggests a halo metallicity [Fe/H] ^^-1.5, in excellent agree- 
ment with measured h alo metallicity distributions, which peak 
at [Fe/H] - -1.6 (e.g. jLaird et al.ll 199811 Venn et al.ll2004l) . 

2.2. Hess diagrams of structure in the stellar halo 

One of the main goals of this paper is to explore the degree 
of substructure in the stellar halo of the Milky Way. One way of 
visualizing this issue is through the inspection of Hess diagrams 
where pairs of lines of sight are subtracted, following Eqn. 
We have done this exercise for three such lines of sight in Fig. 
[T] where we have chosen three line-of-sight pairs where the 

^An extension of this methodology was used bv lXu et alj tlOOl , who use 
the SDSS DR4 to study stellar halo structure using star counts and color distri- 
butions of stars at Galactic latitudes 6 > 55°. 
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Fig. 2. — A histogram of the absolute magnitudes of stars with 0.2 < g—r < 
0.4 in three globular clusters: Pal 5, NGC 5053, and NGC 5024. These dis- 
tributions give an impression of the actual convolution kernel suffered by the 
0.2 < g — r < 0.4 MSTO stars in the halo of the Milky Way when going from 
distance to apparent magnitude. In this work, we choose to approximate this 
distribution for modeling the stellar halo with a Gaussian distribution with 
(Mr) = 4.5 mag and aM,. = 0.9 mag, an appropriate choice for a stellar popula- 
tion with [Fe/H]~ -1.5. 



subtraction should have been close to zero, if the stellar halo 
of the Milky Way were symmetric and smooth. 

The lower left panel of Fig. [T] shows the difference of two 
different lines of sight {l,b) = (300, 70) — a line of sight towards 
the Virgo overdensity and a part of the Sagittarius stream — 
and (l,b) = (60,70); in a symmetric model such a subtraction 
should come out to zero. The grey levels saturate at deviations 
of ±100%. It is clear that the (/,/?) = (300,70) Hne-of sight 
has strong order-of-unity overdensities at MSTOs fainter than 
r > 21, or distances of > 20kpc assuming a MSTO absolute 
magnitude of Mr ~ 4.5. One can see also a weak sub-giant and 
red giant branch feature at g-r ^ 0.5 and 18 < r < 20, again 
indicating distances > 20 kpc. 

The lower middle shows a line of sight towards (l,b) = 
(44, 37) minus the Hess diagram for stars towards (/,/?) = (15,41). 
This subtraction would be expected to come out close to, but 
not exactly, zero. It would be ideal to be able to subtract off 
the 'correct' pairing of {l,b) = (316,37); however, SDSS has 
not mapped that area of sky owing to its southern declination, 
S = -25. The grey scale saturates at ±50%. There are minor ar- 
tifacts in the subtraction; however, one can clearly see an over- 
density of main sequence stars with a MSTO with r ^ 20.5, 
corresponding to a distance of ^ 16 kpc. 

The lower-right panel shows a line of sight towards {l,b) = 
(167,35) — a line of sight towards part of the Low-Latitude 
overdensity — minus that of (l,b) = (193,35). In a symmetric 
halo this subtraction should be identically zero. The grey scale 
saturates at ±50%. There is a weak MSTO overdensity at r 
19 mag, some ~ 7 kpc from the observer. 

While these lines of sight have been selected to show (vary- 



ing degrees) of halo inhomogeneit}Q, they suffice to illustrate 
two key points. First, the halo is far from homogeneous, 
with strong order-of-unity overdensities as well as weaker ^ 
10-20% features. Second, owing to the partial sky coverage 
of the SDSS, it is difficult to visualize and quantitatively ex- 
plore the structure of the Milky Way's stellar halo using CMD 
subtractions. 

2.3. Main sequence turn-off star maps of the stellar halo 

One more intuitive approach to the distribution of stars in 
the stellar halo is to construct maps of the number of MSTO 
stars in different magnitude (therefore, roughly distance) slices. 
We select MSTO stars with foreground extinction-corrected 
0.2 < g-r < 0.4; this color range was selected empirically 
to encompass the most densely-populated bins of color space 
for the halo MSTO stars in Fig.[T] In 32.11 we showed that in 
such a color bin the average absolute magnitude of the MSTO 
stars in that bin were 4.3 and 4.7 respectively for Palomar 5 
([Fe/H] - -1.4) and NGC 5024 ([Fe/H] - -2.1); accordingly, 
we adopt an average MSTO M, = 4.5 in what follows for stars 
in the color bin 0.2 <g-r < 0.4. Such an absolute magnitude is 
in agreement with model CMDs, which suggest M, = 4.7 ± 0.2 
for stars with metallicities [Fe/H] ^ -1.5 ± 0.5. We make the 
assumption that all stars in the stellar halo are 'old' (i.e., ap- 
proximately the same age as the calibrating globular clusters). 

We show 0.5 mag wide bins of r-band magnitude between 
18.5 < r < 22.5, corresponding to heliocentric distances of 
7 < d/kpc < 40. At such heliocentric distances, the vertical 
distance above the Galactic plane is > 5 kpc along all lines of 
sight, or at > 5 scale heig hts following the ~ 900 pc th ick disk 
scale height estimated bv lLarsen & Humphrevsl (l2003h . Thus, 
the dominant contribution to the MSTO maps is from the stellar 
halo. The resulting Lambert azimuthal equal-area polar projec- 
tions, logarithmically-scaled, are shown in Fig.[3fi 

While one loses the ability to probe for population differ- 
ences in the stellar halo because of the broad color bin adopted 
to derive these densities, it is much more straightforward to 
visualize the distribution of halo MSTO stars using this tech- 
nique. From Figs. [T]and|2] one can see that MSTO stars at a 
single distance will show up in multiple distance bins: the bins 
are 0.5 mag wide, and the RMS of a single distance stellar pop- 
ulation is ^ 0.9 mag. This can be seen easily from inspection 
of some of the 'hot pixels' in Fig. [3] corresponding to known 
globular clusters or dwarf galaxies. These features persist from 
map-to-map despite there being a single population at a unique 
distance, giving a visual impression of the covariance between 
the different maps. 



' Although, in fact, we found it impossible to avoid at least low-level inho- 
mogeneity along any pair of lines of sight. 

■*This presentation is similar to that of Fig. 24 of Juric et al. 1 2007), who 
present this kind of analysis for 20 < r < 21, and Newberg et al. (2007), who 
present a similar diagram with slightly more restrictive color cuts for 20 < g < 
21. 



The stellar halo of the Milky Way 
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Fig. 3. — The stellar halo of the Milky Way as seen by SDSS. The grey scale denotes the logarithm of the number density of 0.2 < g-r < 0.4 stars per square 
degree in eight different magnitude (therefore mean distance) slices; a Lambert azimuthal equal-area polar projection is used. The black areas are not covered by the 
SDSS DR5, and reflect the great circle scanning adopted by the SDSS when collecting its imaging data. Apparent 'hot pixels' are stellar overdensities from globular 
clusters and dwarf galaxies. 
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Focusing on the brightest bins, 18.5 < r < 20, correspond- 
ing to heliocentric distances between ^ 7 kpc and ~ 11 kpc, the 
stellar distribution appears rather smooth, with higher density 
towards the Galactic center and Galactic anticenter. In the case 
of the Galactic center, the interpretation is straightforward: one 
is probing lines of sight which pass ^ 5 kpc from the Galactic 
center, and probe the denser inner parts of the stellar halo. In the 
case of the Galactic anticenter, such a structure is not expected 
in a oblate/triaxial halo model, and recalling the < 1 kpc scale 
height of the thick disk cannot be a thick disk; this is the well- 
known Low-Latitude stream (e.g., Newberg et al. 2002; Penar- 
rubia et al. 12005; Momanv et al. 2006). In this visualization, the 
stream appears to be spread out between a few different magni- 
tude bins: at b < 30° some of that spread may be real, but the 
well-defined structure at {l,b) ~ (165,35) has a relatively nar- 
row distance spread (see the Hess diagram residual in the lower 
right-hand panel of Fig. [T] showing a reasonabl^narrow main 
sequence; see also the discussion in Grillmair|[2006bl) . 

Focusing on the more distant bins, 20 < r < 22.5, corre- 
sponding to heliocentric distances between ~ 14 kpc and ^ 
35 kpc, one finds little contribution from the Low-Latitude 
stream. Instead, superimposed on a reasonably smooth back- 
ground is a prominent contribution from large tidal tails from 
the ongoing interacti on of the M ilky Way with the Sagittarius 
dwarf galaxy (see Bel okurov et al. 2006a, for a much more de - 
tailed discussion). As quantified by Belokurov et al.' ('2006a), 
one can discern a distance gradient in the stream, from the clos- 
est populations towards the Galactic anticenter {l,b)^ (200, 20) 
to the most distant populations towards (l,b) ^ (340,50). 

While it is clear from these maps that the stellar halo of the 
Milky Way is not completely smooth, there is a 'smooth' (i.e., 
not obviously structured) component which dominates these 
maps: if there are variations in this component, these must be 
on spatial scales > 10° on the sky (or scales > 1 kpc at the dis- 
tances of interest for this paper). A number of methods could 
be devised to probe halo structures on such scales. In this pa- 
per, we choose to construct models of a smooth stellar halo to 
represent the Milky Way, and to ask about the fraction of stars 
deviating from this smooth global model as a measure of sub- 
structure in the halo. This exercise is the topic of the remainder 
of this paper. 

3. MODELS OF A SMOOTHLY-DISTRIBUTED STELLAR HALO 

The stellar halo of the Milky Way is modeled using an tri- 
axial broken power-law, where we explore oblate and prolate 
distributions as special cases of triaxial. The minor axis of the 
ellipsoid is constrained to be aligne d with the normal to the 
Gal actic disk (this i s is contrast with iNewberg & Yannvll2005l 
and IXu et al.ll200 ^, who allow the minor axis to vary freely). 
There are 7 free parameters: the normalization A (constrained 
directly through requiring that the model and observations have 
the same number of stars in the magnitude and color range con- 
sidered in this paper), the two power laws and aout, the 
break radius rbieak, b/a, c/a, and the Galactocentric longitude 
of the major axis Lmajor- We adopt a grid search, with between 
4 and 10 values in each parameter of interest, yielding typi- 
cally several hundred to several thousand smooth models to test 
against the data. In what follows, we assume a distance to the 
Galactic center of 8 kpc and a M,- = 4.5 for the MSTO stars with 
0.2 < ^-r < 0.4, with a ctm, = 0.9 mag. 

We choose to define the best fit to be the fit for which the 
RMS deviation of the data a around the model is minimized, 
taking account of the expected Poisson counting uncertainty in 



the model, summed over all bins in /, b, and magnitude: 

i'y-) = - V(A-M,)'- - V(m;-m,)' (2) 



a/total = y; ' , (3) 



where D, is the observed number of main-sequence turn-off 
stars in bin /, M, is the exact model expectation of that bin, M- 
is a realization of that model drawn from a Poisson distribution 
with mean M,, and n is the number of pixels. We could have 
chosen instead to define the best fit by minimizing the reduced 



(4) 



where a} is the Poisson uncertainty of the model We have 
chosen not to do so in this case because we are interested in 
quantifying and placing a lower limit on the structure in the 
stellar halo in this paper, not in finding the 'best' fit to the stel- 
lar halo in a sense (we show in ^that the stellar halo model 
with the lowest \^ has a (j/total that is close to, but slightly 
higher than the stellar halo model with lowest cr/total). For the 
purposes of substructure quantification, u/total has two deci- 
sive advantages. Firstly, unlike x^, cr/total is independent of 
pixel scale0 provided that the substructure in the halo is well- 
sampled by the chosen binning scale. Secondly, for cr/total, the 
contribution of Poisson noise to cr has been removed, leaving 
only the contribution of actual halo structure to the variancqj. 
Thus, even though we have adopted a pixel size of 0.5° x 0.5° 
in what follows (corresponding to > lOOpc scales at the dis- 
tances of interest), our results are to first order independent of 
binning scale (because empirically we find that the vast major- 
ity of the variance is contained on ^^kpc scales and greater). 
We defer to a future paper the exercise of understanding and 
interpreting the scale dependence of stellar halo substructure. 
The main uncertainty in the estimated values of cr/total is from 
the major contributions of a few large structures on the sky to 
cr/total, both through influencing the 'best fit' and through their 
direct contribution to the residuals. Later, we attempt to quan- 
tify this uncertainty through exclusion of the most obvious sub- 
structures from consideration before fitting and estimation of 
cr/total. 

The model parameters (including the normalization) give an 
estimate of the total number of stars in the halo. We cal- 
culate the total number of stars contained in the model with 
Galactocentric radius 1 < rcc/kpc < 40. In order to inter- 
pret this value as a mass, it is necessary to convert the number 
of 0.2 < g-r < 0.4 stars into a mass by calculating a mass- 
to-number ratio. We ado pt an empirical approach, following 
INewberg & Yannvl (|2005|) . Given that the Pal 5 MSTO color 
seems to be a good match to the stellar ha lo MSTO color, we 
use the mass of Pal 5 5000 M©; Oden kirchen et al.ll2002l) 
and the number of stars in Pal 5 above background ('--'1069 
stars with 0.2 < g-r < 0.4) to define a mass-to-number ra- 

''The uncertainty in the model is chosen here because we are evaluating the 
likelihood of the dataset being drawn from the model. 

*The quantity ^ (a-) is inversely proportional to n in the presence of in- 
trinsic structure in the dataset, as is the quantity i "^jDj, thus making cr/total 
pixel scale independent. 

^ We have confirmed by rebinning the data by factors of 16 in area that 
cr/total is indeed independent of pixel scale; thus, the dominant contribution to 
the intrinsic structure of the stellar halo must be on linear scales > 400 pc. 
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Fig. 4. — The cr/total of a large number of oblate halo models. Each point 
represents the value of cr/total for a different oblate halo model: open dia- 
monds show the residuals when no clipping is applied, points show the result 
when areas with contributions Sgr/Low-Latitude streamA'irgo are excised be- 
fore carrying out the analysis. In each case, we show the values of cr/total 
as a function of ciout, r^feak and c/a, marginalized over all other model 
parameters. Recall that our definition of cr/total subtracted off the Poisson un- 
certainties already, and is a measure of the degree of substructure on scales 
> 100 pc. It is clear that the oblateness c/a of the halo is the best-constrained 
parameter; combinations of all of the other parameters can provide equally- 
good fits, given an oblateness. Small random offsets are applied to the discrete 
values of «oui, rbrcak and c/a to aid visibility. 



tio - 4.7M0/MSTO staiB This ratio is in excellent agree- 
ment with values derived using stellar population models for 
populations with [Fe/H] < -1.5; these models have values of 
- 4M0/MSTO star 

As is clear from Figs. [3] and [TO] a significant part of the de- 
viations from a smooth stellar halo is driven by the Sagittar- 
ius and Low-Latitude streams, and by the Virgo overdensity. 
We therefore run the whole minimization twice, once allow- 
ing all b > 30° data to define the fit, and a second time mask- 
ing out most of the Sagittarius and Low-Latitude streams, and 
the Virgo overdensity, by masking regions with b < 35° and 
Q < X < 30, where X is the abscissa of the equal-area projec- 
tion: X = 63.63961-\/2(l -sinZ?). This masking is done to con- 
strain the importance of these larger structures in driving the 
model parameters and residual fraction. 

4. RESULTS 

In this section, we present the fitting results for a large set 
of smooth, symmetric stellar halos. In ^4.1| we present the re- 
sults from oblate stellar halos (i.e., the two longest axes have 
equal lengths). In 34.21 we discuss the fitting results for triaxial 
stellar halos (where all three axes can have different lengths), 
comparing this general case to the case of an oblate halo. 

' IKoch et al . 1 2004) find a deficit of low-mass stars in the central parts of Pal 
5, suggesting that this ratio may be a lower limit. 



b 0.025 



— I — ' — I — ' ' — I — ' — r 

Masked Sagiitorius, 
Virgo and Law-Latitude 




— I — ' — I — 
All b>30 



-4 -2 2 4 -4 -2 2 4 
Ap/a hp/o 



Fig. 5. — The distribution of differences between the observed star counts 
per 0.5° X 0.5° pixel and that predicted by the best-fitting smooth model, di- 
vided by the cr predicted by Poisson uncertainties (black lines). The gray line 
shows the expected distribution from Poisson fluctuations around the smooth 
model. The left panel shows the distributions for the case in which sky areas 
of the Sagittarius, Virgo and Low-Latitude overdensities have been excised 
before this analysis; the right panel shows the results for all b > 30° data. 
Note that ~ 1 /2 of the excess variance is in the 'peak' of the histogram (with 
Ap| < 3cr), and the rest of the excess variance reflects a number of pixels with 
Ap| > 3cr (predominantly towards overdensities, rather than towards under- 
densities). 
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Fig. 6. — A comparison between cr/total and for a large number of oblate 
halo models. Each point represents a different oblate halo model: open dia- 
monds show the residuals when no clipping is applied (260456 degrees of free- 
dom), points show the result when areas with contributions Sgr/Low-Latitude 
stream/Virgo are excised before carrying out the analysis (154336 degrees of 
freedom). 
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Fig. 7. — The values of a large number of oblate halo models. The figure 
is formatted identically to Fig.|4] 
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Fig. 8. — Covariance between different model parameters, for the 'best' 
oblate fits (diamonds show models with cr/total < 0.45, whereas crosses show 
models with < 8 X lO') for which all data with b > 30° were fit. Small 
random offsets are applied to the discrete values of ai„, Qout, ^bieak c/a to 
aid visibility. 



In Fig. |4] we show how the residual fraction depends on the 
halo parameters for a survey of parameter space for oblate ha- 
los. It is immediately clear that these smooth models are a 
very ;?oor representation of the structure of the stellar halo, with 
values of tr/total > 0.4 for the best-fitting models for the case 



Fig. 9. — Data fitting results for the triaxial model halos, analogous to Fig. 
|4] In this figure, we show only the behavior of the 'extra' parameters required 
for a triaxial fit, as the behavior of a[^,aaut, ''break ™d c/a is similar to the 
oblate case shown in Fig.|4] Again, diamonds show the results for all data with 
b > 30° , and the points for the case where Sagittarius, the Low-Latitude stream 
and the Virgo Overdensity were masked out. The top two panels show how 
RMS depends on b/a (where b/a = \ is the oblate case and is not shown), and 
^major. the angle between the long axis of and the GC-Sun line. In the bottom 
two panels, we show covariance between L„,iijor and b/a, and b/a and c/a for 
model fits with (T/total< 0.44 (diamonds) and < ^-8 X 10^ (crosses) for 
which all data with b > 30° were been fit. Small random offsets are applied to 
the discrete values of b/a, L,^ajor and b/a to aid visibility. Including triaxiality 
does not significantly improve the quality of fit; when triaxiality is included 
then values of Lmajm between -40 and are favored, reflecting an attempt by 
the triaxial smooth halo model to fit out contributions from the Sagittarius tidal 
stream. 



where all b > 30° data are fit, and a/total > 0.33 for the case 
where Sagittarius, Virgo and the Low-Latitude overdensities are 
clipped. Prolate models were attempted, and were all consider- 
ably poorer fits than the oblate case shown here (i.e., the trend 
towards poorer fits in Fig.|4]with increasing c/a continues for 
c/a > 1). 

In Fig.|5] we show with the black lines the distribution of the 
differences between observed and smooth model distributions 
in 0.5° X 0.5° bins for both the case where Sagittarius, Virgo 
and Low-Latitude structures were masked out (left panel) and 
for all b > 30° data (right panel). In grey, we show the dis- 
tribution expected for Poisson noise around the smooth model 
alone. The difference between the observed histograms and 
the Poisson expectation is the signal which we observe (cr/total 
^ 0.33,0.43 for the clipped and undipped datasets, respec- 
tivelyfl 

From inspection of Fig. |4] it is clear that a variety of differ- 
ent combinations of parameters are able to provide similar val- 

'Note that the appearance of Fig.|5]depends on the adopted binning, through 
the contribution of Poisson uncertainties to the histogram of Ap/cr. The value 
of cr/total is both in principle and in practice independent of binning scale. 
Larger angular bins reduce the contribution of Poisson noise significantly, mak- 
ing the distribution of Ap/cr significantly broader, while the value of cr/total is 
unchanged. 
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ues of cr/total. The oblateness of the halo is best-constrained, 
with values of c/a ^ 0.6 preferrec^B- This determination of 
halo oblatenes s is in excellent agreement with that of previ- 
ous work (e.g.,JChiba & Beers 2000; Chen et al. 20 01; Larsen 
& Humphre ys 2003; Lemon et al. 2004; Newbe rg & Yannvl 
120051 ; IJuric et al.i,2007, ; iXu et aljj20 06). Other parameters are 
less well-constrained: various combinations of a-m, Qfout and 
'"break are Capable of fitting the halo equally well. Best fit stel- 
lar halo masses (over a radius range of 1^0 kpc) come out at 
~ 3.7 ± 1.2 X IO'^Mq for the models with cr/total < 0.45, with 
considerable uncertainty from the mass-to-number ratio. 

In Figs. |6] and|2] we show the relationship between and 
cr/total, and show the run of as a function of the smooth 
halo model parameters. The minimum is 7.65 x 10^ with 
260456 degrees of freedom for the case where all b > 30° data 
are fit, and 3.41 x 10^ with 154336 degrees of freedom for the 
case where Sagittarius, the low-latitude stream and Virgo are 
excised from the fit; in both cases, the probability of the data 
being drawn from the model are zero (to within floating point 
precision). One can see that cr/total and minimizations yield 
similar, but not identical results. The principal difference be- 
tween cr/total and minimization is that models with some- 
what higher c/fl ^ 0.7 are preferred. This is because of the 1 /cr^ 
weighting of y^, that gives higher weight to better-populated 
pixels (in our case, the pixels at larger radii; this tends to give 
Sagittarius high weight in driving the fit). Such a tendency to- 
wards higher c/a with at dist ances >20kpc has been claimed 
before dChiba & Beersll2000l) ; we do not comment further on 
this possible trend here. Nonetheless, the key message of these 
plots is that minimization using x^ and subsequent estimation 
of cr/total yields similar results, but with slightly larger values 
of cr/total than our method, which chooses explicitly to mini- 
mize the metric of interest in order to put a lower limit on its 
value. 

The covariance of the different fitting parameters of the 
oblate case is illustrated in Fig. [8] Models yielding cr/total 
< 0.45 are shown as diamonds, and x^ < 8 x 10^ as crosses, 
where all data with b > 30° are used. It is clear that the degen- 
eracies in Qfin, Qfout and rbieak indicate that there are a number of 
different ways to construct reasonable halo models (see Robin, 
Reyle, & Creze 2000 for similar results), with the general fea- 
tures of a power law aout ~ -3 in the outer parts and a similar 
or shallower power law in the very inner parts of the halo at 
Galactocentric radii rcc ^ 20 kpco- It is important to note that 
the constraints on the 'best fit' halo model are very weak, owing 
to the significant degree of halo substructure. 



angle between the long axis and the line between the Galactic 
Center [GC] and Sun). 

The best triaxial fit is only very marginally better than the 
best oblate fit, with cr/total= 0.42; in particular, the best triaxial 
fit is still a very poor fit to the stellar halo of the Milky Way 
as judged by either cr/total or x^- Inspection of Fig. |9] shows 
that the best models are only mildly triaxial with b/a > 0.8, 
and with Lmajor ~ -20 (roughly lining up with the Sagittar- 
ius stream). In the bottom panels we show the covariance of 
the parameters of all models with cr/total< 0.44 (diamonds) or 
X^ < 7.8 X 10^ (crosses). There is little obvious covariance be- 
tween the 'triaxiality' parameters, or between the power law 
parameters and the triaxiality parameters. This stresses the dif- 
ficulty in fitting a unique model to the halo; owing to the sig- 
nificant degree of halo substructure, there are many ways to fit 
the halo by balancing problems in one part of the halo against a 
better fit elsewhere. 

4.3. A highly structured stellar halo 

The key point of this paper is that a smooth and symmetric 
(either oblate/prolate or triaxial) model is a poor representation 
of the structure of the stellar halo of our Milky Way. The cr/total 
of the b > 30° data around the model is > 42%; even if the 
largest substructures are clipped, the values of cr/total are > 
33% (i.e., the largest substructures contain ^40% of the total 
variance). 

One can obtain a visual impression of how poorly fit the stel- 
lar halo is by a smooth model by examining Fig. [TO] which 
shows the mean stellar density residuals from the best fit oblate 
model. The residuals are smoothed by a 42' Gaussian ker- 
nel to suppress Poisson noise. One can see that the residuals 
are highly structured on a variety of spatial scales. Particu- 
larly prominent are contributions from the well-known Sagit- 
tarius tidal stream (dominating all residuals for 20.5 < r < 
22.5), the Low-Latitude stream (Galactic anticenter direction 
and b < 35°), and the Virgo overdensity (particularly prominent 
in the 19.5 < r < 20 slice as the diffuse overdensity centered at 
(/,/?) -(280,70): se cJuric et al..2007. and ,Newberg et aL2007h . 



4.2. Triaxial models 

The results for triaxial models are shown in Fig. |9l We do 
not show the results for the power-law parameters c^in, c^out and 
'"break, nor the run of cr/total vs. c/a, as the results for these 
parameters is very similar to the oblate halo case. We focus in- 
stead on the results for the 'new' parameters b/a and Lmajor (the 

'"The halo oblateness is alfected by the assumed value of M,-. Variations of 
Mr of ±0.5 mag lead to changes in oblateness of =p0.1. Furthermore, if the stel- 
lar halo has a binary fraction different from that of the globular clusters used to 
calibrate the absolute magnitude and scatter of turn-off stars, the values for ab- 
solute magnitude and scatter would be affe cted at the < 0.3 mag level, leading 
to modest changes in recovered oblateness ILarsen & Humphrevsll2003h . 

" It is interesting in this context that there have been claims of a break in the 
power law of the stellar halo at roc ^ 20 kpc from counts of RR Lyrae stars 
(see iftestoii_e^aLj|1991i although other analyses see no evidence for a break, 
e.g. Jchiba & BeersHlOOO) . 




Fig. 10. — Residuals of the mean stellar density (data— model) from the best oblate model (Qin,Qoui, ''break >c/a) = (— 2.2,-3.5,20kpc,0.7). The panels show 8 
different distance slices, and have been smoothed using a it = 42' Gaussian. The gray scale saturates at ±60% deviation from the model density, and white represents 
an observed excess over the smooth model prediction. 
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There are a number of other less obvious structures. In the 
last three magnitude bins , one can discern the 'Orphan Stream' 
dBelokurov et al.ll2006bHGrillmaidl2006a ). starting at (l,b) 
(250,50) and stretching to {l,b) (170,40) before disappeai-- 
ing into the noise (there is a clear distance gradient, such that as 
I decreases the distance increases). Visible also is a recently- 
identified structure of stars stretching from {l,b) ~ (180,75) 
towards {l,b) ^ (45,45). This structure, called the Hercules- 
Aquila Overdensity by .Belokurov et al.. (2007), extends south 
of the Galactic plane (as shown in that paper) and is at a dis- 
tance of ^ 16kpc from the Sun. The Hercules-Aquila Over- 
density is reflected as a distinct overdensity in color-magnitude 
space, shown in the lower middle panel of Fig. [T] This CMD, 
obtained by subtracting a background field at {l,b) ^ (15,45) 
from an overdensity field at il,b) ~ (44,40), shows a somewhat 
broadened MSTO with turn-off color g-r ^ 0.3 (i.e., a simi- 
lar color to the rest of the stellar halo). Fig. [TO] illustrates that 
this very diffuse overdensity lies in a 'busy' area of the halo, 
making its extent difficult to reliably estimate. There are other 
potential structures visible, in particular in the most distant 
22 < r < 22.5 bin. Some of the structure has low-level strip- 
ing following the great circles along which the SDSS scanri, 
indicating that the structure is an artifact of uneven data qual- 
ity in different stripes. Other structures have geometry more 
suggestive of genuine substructure; we choose to not speculate 
on the reality (or 'distinctness') of these structures at this stage 
owing to the decreasing data quality at these faint limits. 

4.4. Structure as a function of distance 

The visual impression given by Fig. [TO]suggests an increas- 
ing amount of deviation from a smooth halo at larger heliocen- 
tric distances. We quantify this in Fig.[TTl where we show the 
cr/total as a function of apparent magnitude for all stars with 
b > 30° (diamonds). While it is clear that the exact values of 
cr/total will depend somewhat on which smooth model happens 
to fit best, the value of cr/total doubles from distances of ~ 5 kpc 
to ~ 30 kpc. From comparison with the case when Sagittarius, 
the Low-Latitude stream and the Virgo overdensity are removed 
before calculation of the RMS, one can see that much of this in- 
crease in RMS is driven by the few large structures; i.e., much 
of the RMS is contained in a few very well-defined structures 
at large radii. 

5 . COMPARISON WITH EXPECTATIONS FROM A ACDM UNIVERSE 

In this paper, we have attempted to fit smooth models to the 
stellar halo of the Milky Way. Models containing 3.7 ± 1.2 x 
1O**M0 in the radial range 1-40 kpc with power-law density 
distributions p ~ r~^ were favored, although all smooth mod- 
els were a very poor fit to the data. We have found that the 
stellar halo of the Milky Way halo is richly substructured, with 
cr/total > 0.4. The fractional amount of substructure appears to 
increase with radius; this increase is driven primarily by a few 
large structures. 

To put our results into a cosmological context, we compare 
the observations to p redictions for stellar halo structure from 
appropriate models. iBullock. Kratsov. & Weinberg! (12001 1) and 

'^This stri ping has a modest effect on our measurement of cr/total, as illus- 
trated in Fig. llll There are two main effects, working in counteracting direc- 
tions: on one hand, the striping will introduce a small amount of excess vari- 
ance, on the other hand, galaxies misclassilied as stars are smoothly distributed 
across the sky, reducing the variance. We chose to include the 22 < r < 22.5 
bin in the analysis, noting that its exclusion does not affect our results or con- 
clusions. 
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Fig. 11. — The substructure in the Milky Way stellar halo, compared to 
predictions from cosmological models. The cr/total as a function of apparent 
magnitude (distance assuming M,- ^ 4.5) for the 'best fit' oblate model. Dia- 
monds denote the SDSS results for all b > 30° data; crosses denote analogous 
results when the bulk of the Sagittarius and Low-Latitude tidal streams, and 
the Virgo overdensity, have been excised from consideration. The ensemble 
of solid gray lines show the predictions for cr/total from 1 1 models of stellar 
halo growth in a cosmological context from Bullock & Johnston ( 200^; dotted 
hnes are used at small radii where the simulations are hkely to be less robust. 
In these simulations the entire halo arises, by model construction, from the 
disruption of satellite galaxies. 



iBuUock & JohnstonI ( l2005h studied the structure of stellar halos 
created exclusively through the m erg ing and disruption of rea- 
sonably realistic satellite galaxiesQ These studies found that 
the debris from disrupted satellite galaxies produced stellar ha- 
los with: i) roughly power-law profiles with a ^ -3 over 10- 
30 kpc from the galactic c enter (e.g.. Fig. 9 of Bullock & John- 
ston 2005, see also Diemand . Madau. & Moorell2005L Moore 
et al. 2006i) . ii) total stellar halo masses from ~ lO^M© (inte- 
grated over all radii), and Hi) richly substructured halos with 
increasingly evident substructure at larger distances (e.g.. Figs. 
13 and 14 of .Bullock & Johnstoa.2005,) . 



' lAbadi. Navarro. & Steinmetzl 120061) analyzed the properties of the stellar 
halo of a disk galaxy formed in a self-consistent cosmological simulation. Such 
a self-consistent simulation does not require that stellar halos be built up solely 
throu gh accretion; yet, the final halos produced were very similar to those of 
lijullock. Kratsov. & Weinberg. 1 .2001.) and .Bullock & Johnstoa 1 , 2005.) . 




Fig. 12. — Residuals (simulation— smooth model) smoothed using a cr = 42' Gaussian from the best oblate model for Model 2 from lBullock & JohnstonI 120051) in 
8 different distance slices. The gray scale saturates at ±60% from the smooth model density. 
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Fig. 13. — Residuals (SDSS or simulatio ns minus the s mooth model) smoothed using a = 42' Gaussian from the best oblate model fits for the SDSS data (top 
left panel) and for the 11 simulations from lBullock & John ston (2005,). We show only the 20 < r < 20.5 slice, corresponding to heliocentric distances ~ 14kpc. 
The gray scale saturates at ±60% from the smooth model density. 
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Fig. 14. — The fraction of the SDSS footprint in areas with < 70% of the 
smooth model density (upper panel) and > 130% of the smooth model density 
(lower panel) in both the observ ations (blac k line with symbols) and the 1 1 
simulations from iBullock & JohnstonI j2005l gray lines). Dotted lines denote 
where the simulations are argued to be less robust. 



5.1. A quantitative comparison with simulated stellar halos 

We quantify the last statement through comparison of the 
SDS S data for the stellar halo w ith 1 1 simulated stellar halos 
from iBullock & JohnstonI (|2005|H These 1 1 simulated halos 
were generated at random using semi-analytic merger trees ap- 
propriate for a ACDM cosmology for a Milky-Way mass dark 
matter halo. Maps of MSTO stars (analogous to our SDSS data) 
were constructed from the simulated A^-body stellar halos, ac- 
counting for all important observational effects, as follows. The 
number of MSTO stars per particle was estimated using a ratio 
of 1 main sequence turn off star for every 8 Lq, as calibrated 
empirically using Palomar 5. MSTO stars were distributed in 
space by smoothing over the 64 nearest A-body particle neigh- 
bors, using a Epanechnikov kernel of the form (1 -r^). Each 
star was assigned a simulated Galactic latitude, longitude, and 
heliocentric distance assuming that the Sun is 8 kpc from the 
Galactic center. The heliocentric distance is used to generate 

''^The number_of^2ticlesjnJhe_stell^halo^^ Abadi, Navarro, & Stein- 
metz ( 2006) model galaxy was unfortunately too small to permit a proper com- 
parison with the SDSS data. 



r-band apparent magnitudes for the MSTO stars assuming an 
absolute magnitude Mr = 4.5 and scatter (Jm, = 0.9 (following 
@. The models were then placed in a Lambert equal-area pro- 
jection, and the survey limits of the SDSS DR5 data analyzed 
in this paper applied to the simulated maps. These simulations 
were analyzed in the same way as the SDSS data, by fitting the 
same grid of oblate models. The results are shown in Fig. [TTl 
andFigs.[T2]and[T3l 

Fig.[TT]shows the main result of this analysis: all simulations 
predict a great deal of halo substructure, with values of a/total 
> 0.2. The typical smooth halo fitting parameters (where we 
quote the average and scatter derived from the fits to the 11 
simulated stellar halos) are similar to that of the Milky Way's 
halo with Oout ^ -3.4±0.6,Mi<,./kpc<4o 2.8 ± 1.5 x IO'^Mq, 
and c/a^ 0.65 ±0.25; values of ah, within rbieak ^ 25 kpc tend 
to be higher than that observed for the Milky Way at -1 .3 ± 0.7. 
At small Galactocentric radii < 15 kpc, the simulations are ex- 
pected to be much too structu red (owing to the lack of a Uve 
Galactic potential, see §4.2 of iBuUock & Johnstonll200"5h : ac- 
cordingly, we show results for heliocentric distances < 10 kpc 
as dotted lines, and place little weight on the relatively high val- 
ues of ctin recovered by the best-fitting models. At larger radii, 
where the simulation results are expected to be more robust, 
there are model halos with both less structure and more struc- 
ture than the Milky Way's stellar halo. We illustrate this result 
in Figs.[T2]and[T3] Fig.ll2lshows the residuals (simulation-best 
fit smooth halo) for a model with very similar a/total to the 
Milky Way on the same grey scale used for Fig. [10] in eight 
different apparent magnitude slices. Fig. 1 131 illustrates the di- 
versity of simulated halos, showing the 20 < r < 20.5 appar- 
ent magnitude slice (corresponding to heliocentric distances 

14 kpc) for the SDSS and the 1 1 ACDM realizations of Milky 
Way mass stellar halos. A number of the general characteristics 
of the simulations match the characteristics of the SDSS data: 
the angular extent of 'features' in the nearest bins is typically 
very large, whereas the angular width of streams in the distant 
bins tends to be smaller In the distant bins, the halo substruc- 
ture is a combination of well-confined, relatively young streams 
and diffuse sheets of stars from both older disruption events 
and young events on almost radial orbits (K. Johnston et al., in 
preparation), with large-scale overdensities and underdensities 
being seen. 

In Fig. [14] we explore the fraction of area in under- and over- 
densities in both the observations (black lines and symbols) and 
the 1 1 ACDM realizations of Milky Way Mass stellar halos 
(gray lines). We quantify this by comparing the fraction of 
area for each apparent magnitude slice with densities 30% or 
more below the smooth model at that radius (/<o.7, shown in 
the upper panel), and the fraction of the area in each slice with 
densities 30% or more above the smooth model at that radius 
(/>i.3, in the lower panel). This comparison is sub-optimal in 
the sense that both the model and data have a non-zero contri- 
bution from Poisson noise (the immunity to Poisson noise was 
one of the key advantages of the cr/total estimator), although 
we have reduced the Poisson noise by rebinning the data and 
models in 4 x 4 pixel bins; with this rebinning, the variance 
from counting uncertainties is 16 times smaller than the intrin- 
sic variance. One can see the expected result that much of the 
sky area is covered in underdensities, with a smaller fraction of 
the sky in overdensities. Again, the models at Galactocentric 
radii > 15 kpc (where they are reliable) reproduce the general 
behavior of the observed stellar halo rather well. Interestingly, 
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the fraction of sky area in overdensities tends to be somewhat 
lower in the models than in the observations (i.e., there may be 
room for the models to predict niore substructure). 

This comparison shows that the the overall level of the sub- 
structure seen in the Milky Way's stellar halo falls into the mid- 
dle of the range of predictions from simulations — simulations 
in which the stellar halo arises exclusively from the merging 
and disruption of satellite galaxies. Furthermore, the charac- 
ter of the structures in the simulated stellar halos is very simi- 
lar to those observed in the Milky Wa\F^ The models clearly 
have some shortcomings; in particular, the use of a slowly- 
growing rigid potential for the cent ral disk galaxy in th e Bul- 
lock & Johnston (i2005h simulations leads to excess structure in 
the central parts. Furthermore, it is possible that the real stel- 
lar halo has a 'smooth' component formed either in situ in the 
potential well of the galaxy or accreted so early that no spa- 
tial structure remains. Our analysis shows that there is no need 
for such a smooth component to explain the data, and suggests 
that a smooth component does not dominate the halo at radii 
5 < rGc/kpc< 45. Yet, we have not tested quantitatively how 
large a smooth component could lie in this radial range and still 
lead to the observed RMS; such an exercise wiU be the object 
of a future work. 

5.2. Limitations of this comparison 

While there are steps which can and will be taken with this 
dataset to sharpen the comparison with the simulations (e.g., 
a quantitative comparison of the morphology and spatial scale 
of substructure, and the investigation of substructure metallic- 
ities), it is nonetheless clear that 'small number statistics' is a 
key limitation of this work. The SDSS DR5 contiguously cov- 
ers only 1/5 of the sky, encompassing some 5-10% of all halo 
stars, with Galactocentric radii between 5 and 45 kpc (as esti- 
mated by comparison of the smooth halo stellar masses with the 
actual mass contained in the maps). Larger and deeper multi- 
color imaging surveys will be required to expand the coverage 
of the Milky Way's stellar halo, probing to larger halo radii 
where models predict that halo substructure should be easier 
to discern (see, e. g., the prominent substructures discovered by 
ISesar et aTl 120071 using RR Lyrae stars in the multiply-imaged 
'Stripe 82' of the SDSS). Yet, there is significant halo-to-halo 
scatter in the simulated stellar halos; thus, matching the proper- 
ties of a single stellar halo will always be a relatively easy task. 
More powerful constraints will come from studies of the stel- 
lar halos of statistical samples of galaxies using high-resolution 
ground-based o r HST data (see encouraging progress from e.g., 
iFerguson et aLll2002l and lde Jong et al.li2007h . 

6. CONCLUSIONS 

In this paper, we have quantified the degree of (sub-)structure 
in the Milky Way's stellar halo. We have used a sample of 
stellar halo main sequence turn-off stars, isolated using a color 
cut of 0.2 <g-r < 0.4, and fit oblate and triaxial broken power- 
law models of the density distribution to the data. 

We find that the 'best' fit oblateness of the stellar halo is 
0.5 < c/fl < 0.8 over the Galactocentric radial range 5 to 40 kpc. 
Other halo parameters are significantly less well-constrained; 
many different combinations of parameters (including mild tri- 
axiality) can provide comparably good fits. A single power 

'^Model 2 in particular matches tlie trend in RMS witli Heliocentric distance 
and the fraction in over/underdensities to within < 0. 1 for all bins with distance 
> 10 kpc. ^ 



law p (X with a = -3 provides an acceptable fit. Values of 
-2 > a > -4 are also reasonable fits, as are halo profiles with 
somewhat shallower slopes at r < 20 kpc and steeper slopes out- 
side that range. The halo stellar mass of such models between 
Galactocenti-ic radii of 1 and 40 kpc is 3.7± 1.2 x lO^M©, with 
considerable uncertainty from the conversion of the number of 
0.2 < g-r < 0.4 turn-off stars to mass. 

Importantly, we find that all smooth models are very poor fits 
to the spatial distribution of stellar halo stars. Deviations from 
smooth parameterized distributions, quantified using the RMS 
of the data around the model fit in 0.5° x 0.5° bins (>100pc 
scales at the distances of interest) give cr/total > 0.4, after sub- 
tracting the (known) contribution of Poisson counting uncer- 
tainties. Furthermore, the halo seems significantly more struc- 
tured at larger radii than in the inner ^ 10 kpc; a few individual 
structures dominate this increase in cr/total at larger radii. 

Qualitatively, these results show that the stellar 'substructure' 
found in the Milky Way's halo is not at all a small perturbation 
on top of a smooth halo. Remarkably, this same conclusion 
holds when excising the most prominent known substructures 
from the analysis, such as the Sagittarius stream, and then con- 
sidering the remaining area of the sky. 

We compared these observational results with models of stel- 
lar halo growth in a cosmological context taken from Bullock 
& Johnston (2005). In these models, the stellar halo arises 
exclusively from the disruption of and mergers with satellite 
galaxies. The models were analyzed in the same way as the 
SDSS data. Their models predict a ~ -3 in the radial range 
10-30 kpc, halo masses ~ IO^Mq integrated over all radii (or 
masses ~ 3 x lO^M© in the radial range 1-40 kpc), and richly- 
structured stellar halos with cr/total > 0.2. At radii where the 
model predictions are most robust, the models show a range 
of degrees of substructure, from substantially less than that ob- 
served for the Milky Way to substantially more. Furthermore, 
the character of the substructure appears very similar to that 
showed by the Milky Way's stellar halo. While it is clear that 
the models are not perfect, this comparison lends considerable 
quantitative weight to the idea that a dominant fraction of the 
stellar halo of the Milky Way is composed of the accumulated 
debris from the disruption of dwarf galaxies. 
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